
EstimatesResults  = 'Results/full.mat';
SaveCFto = 'Results/full_CF.mat';

load(EstimatesResults);

% Result file must be loaded.

%% Compute long-run of current process:
% Prepare data, states, pproc etc based on options from result file 
% (typically result file will not include this data structure)
optCount = opt; optCount.SelSubSample = 0; 
optCount.TrunkAge = 0; % Option to trucate age (added later, so may not be present in some results files)
[FID_Grid,Action,Age,AgeState,States,pproc,DS] = prepdataM10(optCount);
% Get Baseline value function
[V,Pa,F] = ValueFunction(States,pproc,theta_est,opt);
% Create boxes for sugarcane proximity:
BoxSize = 2000*States.SRDMidPoints(1:end-1,1)';
[Boxes] = BoxesCreate([DS.CoordX,DS.CoordY], BoxSize);

[UHprob,~] = UpdateUHProb2(p_unob_het,States,AgeState,Action,Pa);

IniAgeStates = double(AgeState(:,end)); IniExoTimeStates = double(States.ExoTimeStates(:,end));

Ns = 80;

CF.FID_Grid = FID_Grid;

disp(' Running simulations for baseline case: ')

for s=1:Ns
    [ CF.OutputPathBase(:,s),  CF.YieldPathBase(:,s), CF.AcreagePathBase(:,s), CF.LandUseBase(:,s), CF.OutputDetailBase(:,s)]  = SimulationPathsLong(IniAgeStates, IniExoTimeStates, Pa, States,Boxes, pproc, theta_est, UHprob, opt );
    ProgressBar(s,Ns)
end

CF.LandUseBase = mean(CF.LandUseBase,2);
CF.OutputDetailBase = mean(CF.OutputDetailBase,2);

fprintf('\n');

save(SaveCFto, 'CF', '-v7.3');


%% Price change counterfactual
% Get States and pproc for changed price process:
PriceChange.on = 1;
PriceChange.factor = 1.05;
[StatesChange,pprocChange] = GenStateVariables(DS,optCount,double(AgeState),PriceChange);
% Get price change value function
[VChange,PaChange,FChange] = ValueFunction(StatesChange,pprocChange,theta_est,optCount);

disp(' Running simulations for 1% price change: ')

for s=1:Ns
    [ CF.OutputPathCF05(:,s),  CF.YieldPathCF05(:,s), CF.AcreagePathCF05(:,s), CF.LandUseCF05(:,s), CF.OutputDetailCF05(:,s)]  = SimulationPathsLong(IniAgeStates, IniExoTimeStates, PaChange, StatesChange, Boxes, pprocChange, theta_est, UHprob, opt );
    ProgressBar(s,Ns)
end

CF.LandUseCF05 = mean(CF.LandUseCF05,2);
CF.OutputDetailCF05 = mean(CF.OutputDetailCF05,2);

fprintf('\n');

save(SaveCFto, 'CF', '-v7.3');
